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Abstract 

We study the stability of superfluid Fermi gases in deep optical lattices in the BCS~Bose-Einstein 
condensation (BEC) crossover at zero temperature. Within the tight-binding attractive Hubbard 
model, we calculate the spectrum of the low-energy Anderson-Bogoliubov (AB) mode as well as 
the single-particle excitations in the presence of superfluid flow in order to determine the critical 
velocities. To obtain the spectrum of the AB mode, we calculate the density response function in 
the generalized random-phase approximation applying the Green's function formalism developed 
by Cote and Grifhn to the Hubbard model. We find that the spectrum of the AB mode is separated 
from the particle-hole continuum having the characteristic rotonlike minimum at short wavelength 
due to the strong charge-density- wave fluctuations. The energy of the rotonlike minimum decreases 
with increasing the lattice velocity and it reaches zero at the critical velocity which is smaller than 
the pair breaking velocity. This indicates that the superfluid state is energetically unstable due 
to the spontaneous emission of the short-wavelength rotonlike excitations of the AB mode instead 
due to pair-breaking. We determine the critical velocities as functions of the interaction strength 
across the BCS-BEC crossover regime. 

PACS numbers: 
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I. INTRODUCTION 



The recent realization of superfluidity in Fermi gases [l-6| has opened a new research 
frontier in ultracold atoms 3] • A great experimental advantage of this system is the ability 
in controlling atomic interactions using a Feshbach resonance |8|. This allows us to access the 
crossover between the Bardeen-Cooper-Shrieffer (BCS)-type superfluidity and Bose-Einstein 
condensation (BEC) of bound molecules, which is referred to as the BCS-BEC crossover 



16|. The study of superfluid Fermi gases in the BCS-BEC crossover is expected to offer new 



insights into the phenomena of superfluidity and superconductivity, which can be applied in 
various fields such as condensed-matter physics, nuclear physics, and particle physics. 

One of the most dramatic features of a superfluid system is the dissipationless superfluid 
flow In particular, critical velocities of superfluid flow have attracted much interest in 
various systems such as superfluid '^He [isl, superfluid ^He 19 1, and atomic Bose-Einstein 



condensates 



20 



2l|. It is well known that the underlying mechanisms for the instability of 



dissipationless flow are different in the BCS and BEC regions in a uniform system. Namely, 
the instability of BCS-type superfluids is considered to be dominated by Cooper pair break- 
ing [22|, whereas the instability of Bose superfluids is induced by spontaneous emission of 



phonon excitations [23[ . It is of interest to study how the mechanism of the instability in 
superfluid Fermi gases changes in the BCS-BEC crossover. 

Recently, Miller et al. investigated experimentally the stability of superfluid flow in Fermi 
gases in shallow one-dimensional (ID) optical lattices across the BCS-BEC crossover 24 1. 
They measured superfluid critical velocities, at which the number of condensed atoms starts 
to decrease, by moving the optical lattice potential through the atomic cloud for different 
values of interatomic interaction and lattice depth 2J]. The measured critical velocities 



showed a crossover behavior between the BCS and BEC regimes taking a maximum value 



at the crossover regime 



24l ]. Critical velocities in superfluid Fermi gases in the BCS-BEC 



crossover have been also addressed theoretically in several papers 251-127 L The observed 



25 



261 . However, 



crossover behavior of the critical velocities has been predicted in Refs. 
most of the theoretical papers are limited within a uniform system 25| or a system in 
the presence of a single potential barrier 26|, which cannot be directly compared to the 
experiment using optical lattices in Ref. {24 1. In Ref. j^^], sound propagation in superfluid 
Fermi gases in optical lattices has been studied using the hydrodynamic approximation. 
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However, microscopic calculation of the critical velocities of superfiuid Fermi gases in optical 
lattices has not been worked out yet. 

In this paper, we study the stability and critical velocities of superfiuid Fermi gases in deep 
one-dimensional, two-dimensional (2D), and three-dimensional (3D) optical lattices in the 
BCS-BEC crossover at zero temperature. We a pply the generalized random-phase approxi- 
mation (GRPA) developed by Cote and Griffin |28(] to the attractive tight-binding Hubbard 
model in order to calculate the excitation spectra in the presence of a moving optical lattice. 
For the stability of Fermi gases in the BCS-BEC crossover, two kinds of excitations play 
crucial roles. One is the single-particle excitation which arises when Cooper pairs are bro- 
ken. The other is the collective density-fluctuation mode, the so-called Anderson-Bogoliubov 
(AB) mode 0, [s^. In a uniform system, the single-particle excitation induces the instabil- 
ity of superfluid flow in the BCS regime, while phonon excitation of the Bogoliubov mode, 
which corresponds to the AB mode in the BCS regime, induces the instability in the BEC 
regime 2J-|26|. We flnd that in deep ID, 2D, and 3D optical lattices, the excitation spec- 
trum of the AB mode has a characteristic rotonlike structure and lies below the particle-hole 
continuum due to the strong charge-density-wave (CDW) fluctuation. The energy of the ro- 
tonlike minimum decreases with increasing the superfluid velocity and it reaches zero before 
the particle-hole continuum does, i.e., before pair breaking occurs. As a result, in contrast 
to the uniform case, the instability of superfluid flow in ID, 2D, and 3D optical lattices is 
induced by the rotonlike excitations of the AB mode rather than by pair-breakings. We 
calculate the critical velocities at which spontaneous emission of the rotonlike excitations 
occurs as functions of the interaction strength in the entire BCS-BEC crossover regime. 

This paper is organized as follows. In Sec. [Tll we present the model and formalism. 
We introduce the tight-binding Hubbard model and the Green's function formalism for 
the GRPA. In Sec. IIIIl we present the results for the stability and the critical velocities of 
superfluid Fermi gases in ID, 2D, and 3D optical lattices. We calculate the excitation spectra 
and determine critical velocities as functions of the attractive interaction. We summarize 
our results in Sec. HVl 
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II. MODEL AND FORMALISM 



In this section, we summarize the Green's function formahsm apphed to an attractive 
Hubbard model. This is necessary for the calculation of response functions in GRPA. The 
excitation spectra of collective modes can be obtained as the poles of the response functions. 
Since our major interest is in the AB mode, we calculate the density response function 
assuming an external field coupled with density. To discuss the stability of superfluid states, 
we extend the previous work for the ground state 



31 



321 ] to the current-carrying states. 



A. Green's function formalism 



We consider two-component atomic superfluid Fermi gases with equal populations loaded 
into optical lattices. We suppose that the optical lattice potential is moving with a constant 
velocity —v in the laboratory frame. If the velocity of the lattice potential does not exceed 
the critical velocity, the Fermi gas remains stable in the laboratory frame due to its super- 
fluidity. This situation can be described equivalently in the frame fixed with respect to the 
lattice potential as a superfluid Fermi gas flowing with a constant quasimomentum 2mv, 
where m is the mass of a fermion. In the following, we describe the system in the frame 
fixed with respect to the lattice potential. Namely, we assume a time-independent lattice 
potential and a supercurrent with the quasimomentum 2mv. 

We assume that the optical lattice potential is sufficiently deep so that the tight-binding 
approximation is valid. Thus, the system can be described by a single-band Hubbard model 
as (we set h = k-Q = 1) 

H=-J J2 + H.C.) +UJ2 4t4lCilCit -f^Yl 4aCia , (1) 

{i,j),(J i i,<J 

where Cj^ is the annihilation operator of a fermion on the jth site with pseudospin a =t, i. 
Here, J is the nearest-neighbor hopping energy, U is the on-site interaction energy, and /i is 
the chemical potential. We assume an attractive interaction between atoms {U < 0). 

In order to calculate the density response function, we introduce a fictitious time- 
dependent external field Pj{t) which is coupled with the density. The Hamiltonian with 
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the external field is given by 

K{t) = H + Vit), (2) 

Vit) = 5^P,(t)n„ (3) 
j 

where Uj = cj^c^o- is the number operator. The density response function is obtained by 
taking a functional derivative of the single-particle Green's function by the external field. 
This will be carried out in Sec. Ill Ci 

We use the imaginary time Green's function technique 



The Heisenberg representa- 



tions of the annihilation and creation operators in the imaginary time r are defined as, 

c,,(r) = exp (^^ dr' K{t')^ c.^exp ^ dr' K{t')^ , (4) 

4(r) = exp dr'K{r')^ cj^exp dr' K(r')) . (5) 

We introduce the normal and anomalous single-particle Green's functions, respectively, as 



34| 



G.„.(r,r') = -(T(c,.(r)4(r'))), (6) 

Mry) = -{TMT)c,,{r'))), (7) 

where T(- ■ ■ ) represents the time-ordering operator with respect to r. Using the Nambu 
representation (351] with 

^jir) = ( 'f^\ ) ; ^]iT) = {4^{r),c,,{r)) , (8) 



the single-particle Green's function can be written in the matrix form as 

G,,(r,r') = -(T(vI/,(r)v[/t(r'))) 

G.,,t(r,r') F,,(r,r') 
i^*(r,r') -G,a(r',r) 

We note that in the absence of the external field, the Green's function at equal sites and 
imaginary times Gjj{T, r) is given by 

G,,(r,r) ^ ;im G,,(r,r') (10) 



(9) 



(11) 



where rijcr = Cja^jcr ^^'i = (^ji^^jt P^i^ annihilation operator which is related to the 

wave function of Cooper pairs, as we discuss below. 

From the equations of motion for Cj^ir) and cj^(r), we obtain the equation for the matrix 
Green's function in Eq. as 



d_ 

dr 



J Gij{r, r') + 2J ^ S^osGmjir, t') 
= 5,,,5(r - r') - Ua;{T{n,{T)^,{T)^]{T'))) + P,{T)a^G,,{r, r'), (12) 
where is the Pauli matrix 

^-(l-i)- '''' 

The non-interacting Green's function G^^- (r, r') satisfies 

(-^ + y^^s) G° (r, t') + 2JJ2 5^,l^^GmAr, r) = 5,,5(r - r'). (14) 
From Eqs. (fT2l) and (IH]), the Green's function satisfies the Dyson equation 



G.,(r,rO = G°(r,r') 



, Jo Jo 

l,in 

I -^0 



(15) 



where (3 = 1/T and T is the temperature. In Eq. (fTSll . the self-energy Sij(r, r') is given by 
E.,(r,r') = -t/aaj^ T rfn (r(n,(r)vl/,(r)v[/J(ri))[G,,(ri, r')]-^ (16) 

Here, we introduced the inverse matrix Green's function [Gij{T,T')]^^ , which satisfies 

V ^dn Gu{r,n){Gi,{n,T')r' = 6,J{t - r'). (17) 

Using Eq. f|T7|) . Eq. f|T5l) can be simplified as 

fe(r, tT' = [GUr, tT' - S.,(r, r') - P,(r)a35,,,5(r - r'). (18) 



To ca' 
imation 



cu 
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ate the self-energy in Eq. ( !T6|) . we use the Hartree-Fock-Gor'kov (HFG) approx- 
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(T(r^Kr2)^^(r)^j(rl))) ~ -{ni{r2))G,,{r,r^) + Ga{r,T2)a^Gi,{T2,T{). (19) 

Thus, the self-energy in the HFG approximation is given by 

S.,(r,r') ^ Sf«(r,r') 

= U[{n,{r))a, - <T3G,,(r, T)a^]5,,,5{T - t'). (20) 

In the presence of supercurrent with velocity Cooper pairs are Bose-condensed into 
the state with the center-of-mass quasimomentum q = 2mv. Since the anomalous Green's 



function Fjj^r, r) can be regarded as the wave function of Cooper pairs [3J], it can be written 

as 

= -("^j) = ^exp(2imt; ■ rj), (21) 

where is the superfiuid gap and rj is the location of the jth site. The exponential factor 
on the right-hand side of Eq. fl2T|) describes the supercurrent with quasimomentum 2mv. 

In the presence of the supercurrent, the normal and anomalous Green's functions can be 
written as 

Gij,a{r,r') = exp{imv ■ rij)Gij^a{r,T'), (22) 
F,,(r, r') = exp(2zmi; ■ R,,)Fij{T, r'), (23) 

respectively, where Vij = ri — Vj is the relative coordinate and Rij = (rj + rj)/2 is the 
center-of-mass coordinate of the Cooper pair. Here, Gij^„{T,T') and Fij{T,T') are functions 
of rij. 

To eliminate the phase factors associated with the supercurrent, it is convenient to intro- 
duce an operator \^j(r) and a matrix Green's function Gijij^r') as, 

(c,t(T)exp(— imv ■ r,) \ , ^ 

f ; =7.^.(r), (24) 



(25) 




where the matrix 7^ for the unitary transformation between G and G is given by 

7. = I ' ^ \- (26) 

exp(zmv ■ Vj 

Using Eq. (ITT]) and fl25p . in the absence of the external field, Gjj{T,T) reduces to 

A„/|t;r 

A:/if/| 1 - {„,,) 

In the following sections, we derive equations for G. 



B. Equilibrium Green's function 

In this section, we calculate the equilibrium Green's function in the absence of the external 
field Pj{t) within the HFG approximation introduced in Eq. f|T9|) . Since Gij{T, r') in Eq. fl25|) 
is a function of Vij, we define the Fourier transform of the Green's function as 

Gij{T, t') = ^ exp[ifc ■ - iunir - T')]Gk{iuJn), (28) 

k,uj„ 

where M is the number of lattice sites and w„ = (2n + l)vr//3 is the Fermi Matsubara 
frequency. We note that Gijij, r') cannot be expanded as Eq. ( 128|) because the phase factor 
in Eq. ( 1231) describing the supercurrent depends on the center-of-mass coordinate Rij. 

From Eq. (fT5|) . we obtain the Dyson equation in Fourier space in the absence of the 
external field as 

Gk{iuJn) = Gl,{iujn) + Gl{iujn)t.k{iu:n)Gk{iu:n)- (29) 

Here, Sfe(«a;„,) is the Fourier transform of Sjj(r, r') = 7jSjj(r, t')7*. In Eq. ( 1291) . the unper- 
turbed Green's function G^f^iiujn) is given by 







(30) 

where ^^k = '^J^ui^ — cos k^d) — /i is the kinetic energy, u is the index for spatial dimension, 
and d is the lattice constant. From Eq. (1T9|) . we obtain the self-energy in Eq. ( I29l) as 

^kitUJn) = . 31 

a: 



In deriving Eq. (13T]) . we shifted the chemical potential by the Hartree-Fock energy nU/2, 
where n is the average number of atoms per site. 

By solving Eq. (12 9 p with the self-energy in Eq. ( 13 ip . we obtain the single-particle Green's 
function as 

^u{^uJn) = + — (32) 



where 



A, 



|2 



Vk 



UkVk 



-UkVk u% 



Wfc = 




_^ ik 
^k 


\Vk\^ = 




8k 


UkV*k = 


28k 





(33) 

(34) 

(35) 
(36) 
(37) 



Here, the single-particle excitation energy is given by 

^k — Vk^ £k, 



Et 



(38) 



where 8k = + |A„|2, = {^k+mv + Cfc-m,«)/2, and rjk = {^k+mv - ^k-mv)/^. Equa- 
tion ( l38l) explicitly shows that E^ depend on the superfiuid velocity v. The single-particle 
excitation spectrum in Eq. ( l38i) is shown in Fig. [T] for different superfiuid velocities. In 
Fig. [H the energy gap becomes smaller as increases. When the energy gap reaches zero, 
pair breaking occurs, i.e., v = Vp^,, where is the pair-breaking velocity 36 1. 

We determine the superfiuid gap A„ and the chemical potential fi by solving self- 
consistently the number equation 

2 



n 



M 



(39) 



and the gap equation 



U 



^- = -mE"'^^^[1-2/(^^^)]- 



(40) 
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FIG. 1: Single-particle excitation spectra (upper curves) and E'j^ (lower curves) in 3D optical 
lattices when the superfluid velocity \v\ is (solid line), 0.2/md (dashed line), OA/nid (dash- 
dotted line), and \VpY,\ (dotted line). Here, the superfluid flows along the (7r,7r,7r) direction and 
|i?pb| = 0.628/md. We set U/J = —6.0, n = 0.5, and = ky = kz- 

Equations ( 139|) and ( 140|) are obtained from the diagonal and off-diagonal elements of the 
Green's function in Eq. (l32l) . Here, f{e) = l/[exp(/3£:) + l] is the Fermi distribution function. 
This scheme of solving Eqs. (13^ and f HU]) self-consistently interpolates the weak-coupling 
BCS limit and strong-coupling BEG limit at low temperature when the fluctuation effect 



due to pairs with finite center of mass momenta can be neglected 10| . Throughout the work 
including the calculation of A„ and fi, we assume T = 0. 

In Eq. (40), A„ depends on superfluid velocity v via u^, ffc, and . To explicitly show 
this, we plot A„ in 3D at T = as a function of \v\{< ji'pbl) in Fig. O Note that when 
= at \v\ = \Vpy,\, the superfluid gap does not vanish (A„ 7^ 0), but the superfluid state 
is destabilized due to pair breaking. 
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FIG. 2: Superfluid gap A„ as a function of ji^l in 3D optical lattices, when the superfluid flows 
along the (7r,7r, vr) direction. Here, |A„=o| = 1-92 J and \vpi^\ = 0.628/md. We set U/J = —6.0 and 
n = 0.5. 

C. Response function 



In this section, we calculate the density response function. The density response function 
can be derived by taking a functional derivative of the density by the external field as 



(41) 



Here, we introduce the three-point correlation function Liji{T,T' ,Ti) as 



Lijiir, r',Ti) 



(42) 



where Gij^r, t') = a^Gijir, r'). Using Eq. (H2l) . the density response function can be written 
as 



Xijir,T) = -{T{6n^{T)6nj{T'))) 



(43) 
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where Srii^r) = riiij) — {riiij)) is the density fluctuation operator and 

Lij{T,T')= hm Liij{T,Ti,T'). (44) 

Ti->r+ 

In deriving Eq. (143!) . we used the functional differentiation of the single-particle Green's 
function by the external fleld as 

^-^^^ = (T(^,(r)^t(^.)^^(^^))) (45) 

We note that Eq. ( I45l) is valid within the linear-response regime. 

We derive the equation of motion for the three-point correlation function L. Differenti- 
ating Eq. ffTTj) with respect to the external fleld, we obtain 

+ V [\t2 GUr,r2)^-^^j^asG^j{T2,T'). (46) 
Using the HFG approximation in Eq. f|T9|) . the three-point correlation function satisfles 
Liji{r, T, Ti) = L^ijiir, r, n) + f/ / dra Gim{r, T2)Gmj{r2, r')Xm«(^2, n) 

(T"2,n), (47) 

™ Jo 



where the lowest-order correlation function is given by L°.;(r, r', ri) = Gii{T,Ti)Gij{Ti,T'). 
Thus, the density response function can be obtained by solving Eq. P7j) which is referred 



to as the GRPA equation 28 1. 



To see the diagrammatic structure of Eq. ( 1471) more clearly, it is useful to rewrite Eq. (147|) 
in terms of the irreducible correlation function 

Liji{T, t', Ti) = ^jiir, t', n) - U'S^ / dT2 Gim^r, T2)Lml{T2, ri)Gmj{T2, t). (48) 

Using Eq. (iHD, Eq. (gTl) reduces to 

Liji{T, r', Ti) = Liji{T, t', ri) + Uy2 / dT2 Lij^^r, r', T2)xmi{r2, n). (49) 

™ Jo 



It is clear from Eq. (148!) that L includes the ladder diagrams. On the other hand, Eq. 
includes the bubble diagrams which lead to the random-phase approximation (RPA) 



9 
28| 
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In a homogeneous system, the contribution from the bubble diagrams for an attractive 



interaction can be neglected 



28 



However, in our lattice system, since the bubble diagrams 
induce the instability due to the CDW fluctuation, it is crucial for the analysis of the stability 
of the system to keep the bubble diagrams in Eq. fH^ . We compare the excitation spectra 
with and without the contribution from the bubble diagrams and give a detailed discussion 
of the effects of the CDW fluctuations on the stability of the system in Sec. IIIII 



We solve Eqs. pHj) and f H9|) to calculate the density response function in Eq. f HT|) . We 



define the Fourier transform of Li At, t') as 



Lij{r, r') = ^ exp[iq ■ - iQ^ir - T')]Lg{iQr^ 



(50) 



where fi„ = 2mi / [5 is the Bose Matsubara frequency. From Eq. ( H3l) . the Fourier component 
of the density response function is given by x<j(^^n) = L^qii^n) + L'^qii^n)- 
The Fourier transforms of Eqs. f HHj) and f H9|) are represented as 

U 



(3M 



and 



respectively, where 



LqiiVLn) = LqiiVLri) + U Lq{iQn)Xqi'>'^n) 



Lj(iOn) = ^ Gk{iUJn)Gk-qiiUn - i^n)- 



(51) 



(52) 



(53) 



k.LOn 



In order to rewrite Eqs. (15T|) and (152|) in simpler forms, we define a column vector Cq{iVLn) 
as 



Cq{iVln) 



(54) 



Here, we have used the notation L^^ = L^^iiVLn)- The same notation is adapted for C,q{iVLn) 
and £g(if2n)- In addition, we define a 4 x 4 matrix V as 

/ £)1111 /)1121 /)1211 £)1221 \ 

£)1112 £)1122 £)1212 £)1222 

£)2111 £)2121 £)2211 £)2221 

Y £)2112 £)2122 £)2212 £)2222 y 



Vq{iflr 



(55) 
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where 

^r'(^^«) ^jmT. + (zc.). (56) 

k.UJn 

From Eq. (l56l) . D describes a single bubble diagram of a particle- hole excitation. Thus, 
Eqs. (15T|) and (152|) can be written in the matrix forms as 

= cl{in^)-uv^{inn)c^{inr,), (57) 

Cq{ifln) = Cq{iQn) + UCq{iQn)Xq{i^n), (58) 

respectively. Solving Eq. flFTI) . we obtain 

= [/ + f/:D^(za„)]-i/:°(zfi„). (59) 

Here, / is the 4x4 unit matrix. After the analytic continuation iQn uj + iS (we take the 
limit 6 — )■ +0 after the calculation), we obtain the density response function as 



The excitation spectrum of the AB mode is obtained from the pole of Eq. (1601) . 
If we only take into account the ladder diagrams, Xq(w) reduces to 

Xj(cc;) = [4(a;)]i + [Cqiuj)],. (61) 

Equation (16T1) shows that the pole of x^i'^) coincides with that of Cq{uj) which is obtained 
from the condition |/ + U'Dq{uj) \ = in Eq. fl59|l . 

We calculate the single bubble diagram Vq^u) in Eq. (156|) . Substituting Eq. (l32l) into 
Eq. fl56|) . we obtain 



k 

+^k-q^k ... ^k+q^k , p+ . p+ . 



(62) 



where A = a^A and i? = a^B. In deriving Eq. f l62p . we used EZ^. = ^E^. We assume 
E^ > because we are interested in the stability of Fermi gases before the pair breaking 
sets in, i.e., \v\ < l^pbl- At T = 0, the first and second terms in Eq. (l62l) vanish from this 
condition. The density response functions in Eqs. ( 160|) and (!6T|) are calculated by using 
Eqs. dSnD and dSg). 
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From Eq. ( 162|) . the spectrum of the particle-hole excitation is given by 

<(fc) = K + E\^^. (63) 

For fixed q, uj^{k) makes a continuum for different fc, as shown in Fig. |3l The up- 
per and lower boundaries of the particle-hole continuum are given by miufc and 
maxfc [(jjP'^(fc)], respectively. 



III. RESULTS 



In this section, by calculating the excitation spectra of the AB mode and the single- 
particle excitation, we discuss the stability of superfluid flow, and determine the critical 
velocities of superfluid Fermi gases in ID, 2D, and 3D optical lattices. For this purpose, we 
calculate the dynamic structure factor Sq{u) = —lm[xq{uj)]/iT, which describes the response 
of the system to density perturbations with momentum q and frequency u. The dynamic 



structure factor can be directly measured in experiments by using Bragg spectroscopy 



371, 



Since the GRPA used in this paper is based on a mean-fleld approximation, it is more 
reliable for higher dimensions. Nevertheless, calculations of the excitation spectra in the 
simplest situation of ID can be useful for understanding the essence of the physics governing 
the critical velocity of superfluid fermions in a lattice. Hence, we flrst discuss the excitation 
spectra and the stability of superfluid flow in ID lattices. We note that mean-fleld theories 
have been widely used to qualitatively describe excitations of trapped atomic gases even in 
ID because the flnite size (~ lOOd) speciflc to cold atom systems excludes long- wavelength 



phase fluctuations that destroy the long-range superfluid order |39|-|41|. 

In Fig. [3], we show the dynamic structure factor Sq{uj) in ID optical lattices to illus- 
trate the basic properties of the excitation spectra. One clearly sees that the AB mode 
spectrum lies below the particle-hole continuum. In addition, the AB mode spectrum has 
a characteristic structure with local minima at short wavelengths which is similar to the 



roton spectrum in superfluid ^He [421]. Then, it is expected that as the superfluid velocity 
increases, the energy of one of the rotonlike minima decreases and it reaches zero before the 
lower boundary of the particle-hole continuum does. This indicates that the instability may 
be induced by the rotonlike excitations of the AB mode rather than by the single-particle 
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FIG. 3: Dynamic structure factor Sq{uj) in ID optical lattices. The upper gray region and the 
lower curve correspond to the particle-hole continuum and the AB mode spectrum, respectively. 
The dashed line is the lower boundary of the single-particle continuum. We set n = 0.5 (quarter 
filling), U = — 2.0J, and u = 0. The superfluid gap and chemical potential are |Aj;| = 0.409J and 
^ = 0.624J. Here, in numerical calculation, 5 is set to be small but finite (1.0 x 10"^) so that the 
peak of the AB mode spectrum has a small finite width. 

excitations because the single-particle excitations start to have negative energies when the 



particle-hole continuum reaches zero energy 



25| . Indeed, we will show that this is the case 



in all of ID, 2D, and 3D optical lattices in the remainder of this section. 

We note that this AB mode-induced instability does not occur in superfluid Fermi gases in 
uniform 3D systems. It was found that the instability of superfluid Fermi gases in uniform 
3D systems is induced by pair breaking 
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25| which is associated with the appearance 



of single-particle excitations with negative energies. In this case, the AB mode spectrum 
merges into the particle-hole continuum in contrast with the behavior in Fig. |3] where the 
AB mode spectrum is separated from the particle-hole continuum. As a result, the particle- 
hole continuum reaches zero energy before the phonon part of the AB mode spectrum starts 
to have negative energy as the superfluid velocity increases. This leads to the instability 
induced by single-particle excitations, i.e., by pair breaking. The rotonlike structure of the 
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AB mode spectrum also appears in a uniform ID system 43|] and a 2D lattice system 441.145 1. 

In the following, we first calculate the pair-breaking velocity t'pb- We next discuss the 
behavior of the AB mode spectrum in ID, 2D, and 3D optical lattices by calculating the 
dynamic structure factor Sg{uj). To discuss the excitation spectra, we only show the lower 
boundary of the particle-hole continuum and the peak of the AB mode spectrum in Sq{uj) 
because other details are not necessary for determining the critical velocities. 



A. Pair-breaking velocity 



In this section, we calculate the pair-breaking velocity Vp\^, which can be analytically 



obtained from the condition that the 



zero energy, i.e., minklcu^ (k)] = |25| 



ower boundary of the particle-hole continuum reaches 



In ID case, the pair-breaking velocity is given by 



Vph 



I A, 



sm 



(64) 



md \^ ^^(4J - /i) 

In 2D case, when the supercurrent is flowing in the (vr, vr) and (vr, 0) directions, the pair 
breaking velocities are calculated as 



\v 



pbl 



V2 
md 



sm 



and 



1 • - J 

sm ' 



md 
1 




(65) 



\ V^/i(4J- /i) 
A, 



sm 



(66) 



(/i > 2J), 



md \ 2J 

respectively. We address the stability of superfluid states in these two cases in Sec. IIII CI 
In 3D case, when the supercurrent is flowing in the (tt, tt, tt) direction, the pair-breaking 
velocities are calculated as 



lA, 



|Vpb| = ^sin-^ , ' . (67) 

Since the order parameter and chemical potential depend on the superfluid velocity v, we 
must determine Vp^ by solving Eq. flM|) . (1^ . (1^ . or f E7|) self-consistently with Eqs. (15^ 
and (HUI) . In the BCS limit (|f/| ^ J), t>pb approaches zero because |A„| becomes small. 
On the other hand, if \U\ is so large that |A„| is larger than the denominator in sin~^ in 
Eqs. ( 16^ -( 167|) . Wpb is not definable. 
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FIG. 4: Excitation spectra in ID optical lattices for (a) current-free {v = 0) and (b) current- 
carrying (v = 0.21/?7id) cases. Solid line and dashed line represent the spectrum of the AB mode 
which corresponds to the (5- function peak of Sq{uj) and the lower boundary of the particle-hole 
continuum, respectively. We set n = 0.5 (quarter filling) and U = —2.0J. The superfluid gap and 
chemical potential are (a) \Ay\ = 0.409 J and /i = 0.624 J, and (b) \Ay \ = 0.420 J and n = 0.655 J. 

B. Stability in ID optical lattices 

In this section, we study the stability of superfluid Fermi gases in ID optical lattices. 
In Fig. m we show the dynamic structure factor Sq{uj) in ID optical lattices. It is clearly 
seen in Fig. H] that the AB mode has a gapless and linear dispersion in the long- wavelength 



limit ( 

mode 



qd \ <C 1), which is consistent with the fact that the AB mode is a Nambu-Goldstone 



35|. We obtain the analytic form of the phonon-like dispersion relation in the long- 



wavelength limit {\qd\ ^ 1), as 



d{2J — fi) tan{mvd)q 
1 



a 



a) ta.n^ {mvd)]}\q\ 



(68) 



where vp = 2Jdsm{kpd) cos^mvd), k-p = \ cos ^{{2J — fi)[2J cos{mvd)] ^}\/d, a = NqU, and 
A'o = d{7ivp)^^. The details of the derivation of Eq. ( 168|) are summarized in the Appendix. 
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When V = 0, Eq. (168|) in the BCS hmit (|A^| ^ /i) reduces to the well-known form ojg = 
vf\/1 + NoUlql, which was first obtained by Anderson for a uniform 3D system 29|, |3l| (the 
dispersion has the coefficient l/y/S in 3D case). We note that Eq. ( 1681) is different from 
the dispersion relation obtained in Ref. by using the hydrodynamic and tight-binding 
approximations. 

As we pointed out earlier, in Figs. Hl^a) andlU^b), it is clearly seen that the excitation 
spectrum of the AB mode lies below the particle-hole continuum, and the AB mode spectrum 
has roton-like minima at |g| ~ 2kp [kp = n7T/{2d) is the Fermi wave number in a non- 
interacting ID system]. As v increases, the whole spectrum leans toward the left side and 
the energy of the rotonlike minimum with g < decreases as shown in Fig. Ht^b). As a result, 
at a certain velocity Vc, the rotonlike minimum reaches zero energy, but this occurs before 
';he lower boundary of the particle-hole continuum does. According to the Landau criterion 



23| . this indicates that the spontaneous emission of roton-like excitations of the AB mode 



is induced when v > Vc before pair-breaking occurs at fpb. Thus, the critical velocity is 
given by Vc at which the superfluid flow is destabilized due to the spontaneous emission 
of rotonlike excitations of the AB mode. We note that the phonon part of the AB mode 
spectrum becomes negative at a certain velocity larger than Wpb. This means that phonon 
excitations of the AB mode are irrelevant to the critical velocities. 

This instability driven by negative-energy rotonlike excitations of the AB mode corre- 
sponds to the energetic instability called Landau instability. Another type of instability 
of superfluid called dynamical instability was also proposed for Fermi gases in optical lat- 
tices 271, |46j. Dynamical instability is associated with the appearance of complex energy 
excitations which was flrst observed in Bose condensates in optical lattices 47|. One can 
distinguish the dynamical instability from the Landau instability by identifying the quasi- 
momentum q of the excitations causing the instability. If the instability is caused by the 
excitations at the boundary of the flrst Brillouin zone, e.g., qx = ivr/fi or Qy = ±iT/d in 2D, 
it is the dynamical instability because these excitations inevitably couple with their anti- 



phonon branches 48|, |49|. Indeed, we will see that in 2D lattices the dynamical instability 



due to the AB mode ai q = (ivr/rf, ivr/rf) can occur near the half fllling or in the BEG 
region. 

In Fig. |5l we show Vc and Upb as functions of U/J in the BCS-BEC crossover region. 
One clearly sees that Vc is smaller than Vp^^. We conflrmed that Vc is smaller than Vp^ in 
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FIG. 5: Critical velocity (solid line) and pair-breaking velocity i)pb (dashed line) as functions of 
U/ J in ID lattices. We set n = 0.5 (quarter filling). From Eq. ()64p . fpt, approaches Tr/lmd when 
\U\/J —7- 4 (filled circle). When \U\/ J > 4, the pair-breaking velocity is not definable. 

the entire BCS region (—1 ^ U/J < 0). Thus, the instabihty is always induced by the 
roton-like excitations in these regions. The difference between Vc and Wpb increases with 
increasing \U\/J when one approaches the BEC regime. In addition, both Vc and fpb grow 
monotonically with increasing the interaction \U\/J. 

When \U\/J ^ 1 (BEC region), the size of the Cooper pairs becomes smaller than the 
lattice spacing and each Cooper pair forms a tightly-bound molecular boson. In this region, 
the Hubbard model of Eq. ([T]) is mapped onto a hardcore Bose- Hubbard model with nearest- 
neighbor repulsive interactions (or equivalently the spin-i model) 29|, |4J, |46|, |50| . Since 
it is well known that any kinds of mean-field theory completely fail to describe hardcore 
bosons in ID, our GRPA is also invalid in the BEC limit in ID. Hence, we postpone the 
discussion of the BEC limit to the next section, where we will show results in 2D. 

To discuss the origin of the roton-like minima of the AB spectrum, we show the dynamic 
structure factor when n = 0.9 in Fig. El Compared to the AB mode spectrum when n = 0.5 
in Fig. m the roton-like minima have lower energies than those in Fig. |H It turns out that as 
one approaches half filling (n = 1), the energy of the rotonlike minima becomes smaller. As 
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qdln qd/n 

FIG. 6: Excitation spectra in ID optical lattices for (a) current-free {v = 0) and (b) current- 
carrying (v = 0.12/md) cases. Solid line, dash-dotted line, and dashed line represent the spectrum 
of the AB mode obtained from Eq. (j60p . Eq. (|6ip . and the lower boundary of the single-particle 
excitation continuum, respectively. We set n = 0.9 and U = — 2.0J. The superfluid gap and 
chemical potential are calculated as (a) |A^| = 0.345 J and fi = 1.69 J, and (b) \Ai,\ = 0.350 J and 
/i = 1.70J. 

is well known, the fluctuation due to the formation of charge- density- wave (CDW) order is 



enhanced near half filling in lattice fermion systems [5l|. Thus, the CDW fluctuation leads 
to the rotonlike minima in the AB mode spectrum. At half filling, the rotonlike minima 
reach zero energy even in the current-free case {v = 0) and the superfluid ground state 
becomes unstable due to the formation of CDW order. We show the critical velocity Vc 
when n = 0.9 as a function of f// J in Fig. [71 The critical velocity when n = 0.9 in Fig. [7] 
is smaller than that when n = 0.5 in Fig. [5] due to the strong CDW fluctuation. We notice 
that Vc becomes almost constant below a certain value of interaction {U/J^ —2), which 
reflects the fact that the energy difference between the lower boundary of the particle-hole 
continuum and the roton-like minimum becomes large as \U\ increases. This indicates that 
the CDW fluctuation is enhanced below this value of interaction. 
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FIG. 7: Critical velocity (solid line) and pair-breaking velocity Vp\^ (dashed line) as functions of 
?7/J in ID lattices near half filling. We set n = 0.9. 



To support the above consideration on the origin of the roton-hke minimum, we compare 
the AB mode spectra calculated by Eqs. f lBDl) and ([HID in Fig. O As discussed in Sec. Ill CI 
Eq. (16 ip includes only the contribution from the ladder diagrams, while Eq. (16 Op includes the 
contributions both from ladder and RPA-type bubble diagrams. In Fig. [6l one clearly sees 
that the AB mode spectrum calculated by Eq. fl60l) lies below the one calculated by Eq. f l6T|) 
which actually lies close to the lower boundary of the particle-hole continuum. Since the 
RPA-type bubble diagrams include the effect of the CDW fluctuation [5l|, this behavior of 
the AB mode spectrum is consistent with the above consideration that the CDW fluctuation 
leads to the roton-like structure of the AB mode spectrum. 



C. Stability in 2D optical lattices 

In this section, we discuss the stability of Fermi gases in 2D optical lattices. Here, we 
restrict ourselves to two characteristic situations where the superfluid flows along the (vr, tt) 
or (tt, 0) directions (see Fig.|8]) in order to see the effects of CDW fluctuations on the stability 
of superfluid Fermi gases. 
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FIG. 8: Schematic picture of 2D optical lattices. Arrows indicate the superfluid velocity in (a) 
(vr, vr) and (b) (vr, 0) directions. Each circle represents a lattice site. 



First, we discuss the case when the superfluid flows along the (vr, vr) direction. We assume 
the superfluid velocity v = {v,v)/\/2, where v = \v\. We calculate the dynamic structure 
factor Sq{uj) only when = Qy because superfluid flow is expected to be most unstable for 
excitations with momenta in the opposite direction to the flow. 

In Figs. M^a) andlHl^b), we show the excitation spectra when n = 0.5 and = Qy When 
V = [see Fig. [9]^a)], the rotonlike structure is slightly seen in the AB mode spectrum. As v 
increases, the rotonlike structure becomes remarkable and one of the rotonlike minima goes 
down. At a certain velocity Vc smaller than the pair-breaking velocity fpb, the energy of 
the rotonlike minimum reaches zero [see Fig. [9](b)]. Thus, as ID case, the critical velocity is 
given by Vc at which the instability due to spontaneous emission of rotonlike excitations of 
the AB mode sets in. 



In Figs[9t^c) and[9l^d), we show the excitation spectra in the BEC region [U = — 12J) 52|. 
There we see that the rotonlike minima of the AB mode are present also in this region. 
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FIG. 9: Excitation spectra in 2D optical lattices when the superfluid flows along the (vr, vr) direction. 
Solid and dashed lines represent the spectrum of the AB mode and the lower boundary of the 
particle-hole continuum, respectively. The superfluid velocity is [(a),(c)] u = 0, (b) v = 0.467/?n(i, 
and (d) v = 0.6502/md We set n = 0.5, = Qy, [(a),(b)]C/ = -4.5J, and [(c),(d)]C/ = -12.0J. 
The superfluid gap and chemical potential are calculated as (a) |A„| = 1.33J and /i = 2.23J, 
(b) |A„| = 1.39J and fi = 2.29J, (c) |A^| = 4.92J and = 0.687J, and (d) |A„| = 4.97J and 
/i = 0.746J. 
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FIG. 10: Critical velocity Vc (solid line) and pair-breaking velocity fpb (dashed line) as functions 
of U/J in 2D lattices when the superfluid flows along the (vTjTt) direction. We set (a) n = 0.5 
and (b) n = 0.8. From Eq. Vpb approaches y/2T:/2md when \U\/J — )• 8 (filled circle). When 

\U\/J > 8, the pair-breaking velocity is not definable. 

The critical velocity in the BEC region is also determined by the rotonlike excitations [see 
Fig. [^d)]. Since the roton-like minima are shifted to Qx = Qy = ivr/rf, the instability caused 
by the roton-like excitations is the dynamical instability jSSj . The shift of the roton-like min- 
ima to the edge of the Brillouin zone can be understood as follows. As mentioned before, 
in the BEC region the Hubbard model can be reduced to a hardcore Bose-Hubbard model 
with nearest-neighbor repulsive interactions. The nearest-neighbor repulsion enhances den- 
sity wave fluctuations with the wave vector k = {ji /d^Ti /d) 5J] , leading to the rotonlike 
minimum at q = {tx / d^Ti / d). Thus, the shift of the roton minima means that as one ap- 
proaches the BEC region, the origin of the roton minima changes from the nesting effect 
of the Fermi surface to the nearest-neighbor interactions between molecular bosons. Notice 



that in the 
any longer 



imit of the low filling (n — )■ 0), the roton minima of the AB mode do not survive 



44j and the critical velocity is determined by the long-wavelength part (phonon 



branch) of the AB mode. 
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FIG. 11: Excitation spectra in 2D optical lattices when the superfluid flows along the (vr,7r) 
direction. Solid and dashed lines represent the spectrum of the AB mode and the lower boundary 
of the particle-hole continuum, respectively. The superfluid velocities are (a) f = and (b) 
V = 0.260/?7T,(i. The filling is n = 0.8. We set U = — 4.5J and = Qy The superfiuid gap and 
chemical potential are calculated as (a) |A„| = 1.60J and /x = 3.32J, and (b) |A^| = 1.62J and 
fi = 3.33J. 

In Fig. [TOl we show the critical velocity Vc and the pair-breaking velocity t'pb as functions 
oi U/J. Vc and fpb show qualitatively the same behavior as in ID case. Namely, Vc is smaller 
than Vph and they increase monotonically with increasing the interaction strength \U\/J. 

Near half filling, the rotonlike structure when v = becomes more remarkable due to 
strong CDW fluctuation [see Fig JTTT a)]. Since it is seen in Fig. [TTT b) that the instability is 
driven by the rotonlike excitations with q = (±7r/d, ivr/d), it is the dynamical instability. 
As a result of the strong CDW fluctuation, the critical velocity near half filling is smaller 
than the one at quarter filling, as shown in Fig. [TOT b). As in ID case, the CDW fluctuation 
is strongly enhanced below a certain value of interaction {U/J ^ —3) in Fig. [TOT b). As a 
result, Vc is almost constant below this value of interaction. 

Next, we discuss the case when the supercurrent flows along the (vr, 0) direction in 2D 
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FIG. 12: Excitation spectra in 2D optical lattices when the superfluid flows along the (vr, 0) direction 
for (a) current-free {v = 0) and (b) current-carrying (v = 0.585 /md) cases. Solid and dashed lines 
represent the spectrum of the AB mode and the lower boundary of the particle-hole continuum, 
respectively. We set n = 0.5, U = — 4.5J, and qy = 0. The superfluid gap and chemical potential 
are calculated as (a) |A„| = 1.33J and fi = 2.23J, and (b) |A„| = 1.41J and = 2.31J. 

optical lattices (see Fig. [8]). We assume the superfluid velocity as v = {v,0). The excitation 
spectra when n = 0.5 is shown in Fig. [121 Here, we assume Qy = from the same reason for 
the (vr, vr) case. The behavior of the energy spectra is qualitatively the same as the previous 
cases, i.e., the AB mode spectrum has the rotonlike structure and lies below the particle- hole 
continuum. As the superfluid velocity increases, the AB mode spectrum is pushed down and 
the energy of the rotonlike minimum decreases. The instability sets in at Vc when the energy 
of the roton-like minimum becomes zero. 

In Fig. [131 W6 show the critical velocity Vc and pair-breaking velocity as functions 
of U/J. They also show qualitatively the same behavior as the previous cases. Comparing 
Figs. [T0land[T3l we find that Vc in the (vr, vr) case is smaller than the one in the (vr, 0) case. 
This is because the nesting effect of the Fermi surface occurs in the (vr, vr) direction so that 
it is enhanced by the supercurrent in parallel to this direction. 
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FIG. 13: Critical velocity (solid line) and pair-breaking velocity Upb (dashed line) as functions 
of C//J in 2D lattices when the superfluid flows along the (vTjO) direction. We set n = 0.5. 

We note that near the half filling, the AB mode spectrum with ^j, 7^ [i.e., not in the (tt, O) 
direction] may reach zero before that in the (vr, 0) direction does, even when the superfluid 
flows in the (tt, 0) direction. This is also due to the strong nesting effect in the (vr, tt) 
direction. However, even in this case, our main conclusion remains unchanged. Namely, the 
instability is induced by the AB mode excitations. This effect was also pointed out in the 
BEC regime Q- 



D. Stability in 3D optical lattices 

Let us finally discuss the stability of Fermi gases in 3D optical lattices. We restrict 
ourselves to the situation where the superfluid flows along the (vr, tt, tt) direction. We show 
the excitation spectra in Fig. [T3] for g^. = % = and the critical velocity and the pair- 
breaking velocity as functions of U / J in Fig. [151 It is clearly seen that the AB mode 
lies well below the single-particle continuum and that in the entire region of the attractive 
interaction, the roton part of the AB mode reaches zero before the single-particle continuum 
does. This leads to the conclusion that the critical velocity of superfluid Fermi gases in a 
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deep lattice is determined by the roton part of the AB mode except in the low filling limit, 
regardless of the dimensionality of the system. 



IV. CONCLUSION 

In conclusion, we have studied the stability of superfiuid Fermi gases in ID, 2D, and 3D 
optical lattices at T = 0. By applying the GRPA Green's function formalism developed 
by Cote and Griffin 281] to the attractive Hubbard model, we calculated the excitation 



spectra of the AB mode as well as the single-particle excitation in the presence of superfiuid 
fiow. We found that the AB mode spectrum has the characteristic rotonlike structure being 
separated from the particle-hole continuum due to the strong CDW fiuctuation. The energy 
of the rotonlike minimum decreases as the superfiuid velocity increases and it reaches zero 
at the critical velocity before pair-breaking occurs. This indicates that the instability of 
superfiuid fiow in ID, 2D, and 3D optical lattices is induced by the spontaneous emission of 
the rotonlike excitations of the AB mode. We calculated the critical velocity as functions 
oiU / J and confirmed that it is smaller than the pair-breaking velocity Wpb in the BCS and 
BCS-BEC crossover region. We also found that the CDW instability is strongly enhanced 
near half filling which leads to the suppression of the critical velocity when the attractive 
interaction is large. 

Finally, we remark that our results are valid for superfiuid Fermi gases in deep optical 
lattices because we employed the tight-binding Hubbard model. From this reason, our results 
for ID optical lattices cannot be directly compared to the experiment with shallow optical 
lattices in Ref. 2J|. However, superfiuid Fermi gases have been already achieved in deep 



optical lattices in Ref. 



55| . Our theoretical predictions in this paper may be verified if a 



superfiuid Fermi gas is prepared in a moving deep optical lattice. 
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FIG. 14: Excitation spectra in 3D optical lattices when the superfluid flows along the (vr, vr, vr) 
direction for [(a),(c)] f = 0, (b) = 0.5969/m(i, and (d) v = 0.7894/m(i. Solid and dashed lines 
represent the spectrum of the AB mode and the lower boundary of the particle-hole continuum, 
respectively. We set n = 0.5, Qx = Qy = Qz, [(a))(b)] U = — 6.0J, and [(c), (d)] U = — 14.0J. 
The superfluid gap and chemical potential are calculated as (a) |A„| = 1.92J and fi = 3.78J, 
(b) |A^| = 1.98J and fi = 3.85J, (c) |A„| = 5.70J and /i = 2.11J, and (d) |A^| = 5.77J and 
H = 2.18J. 
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FIG. 15: Critical velocity Vc (solid line) and pair-breaking velocity Upb (dashed line) as functions 
of U / J in 3D lattices when the superfluid flows along the (vr, vr, vr) direction. We set n = 0.5. 
From Eq. ([671), Vpb approaches ^/^^T/2md when \U\/J 12 (filled circle). When \U\/J > 12, the 
pair-breaking velocity is not definable. 

Appendix A: AB phonon spectrum in ID 



Here, we give a detailed derivation of the AB mode spectrum in ID in the long-wavelength 
Umit in Eq. ( I68l) . The matrix T>q[iQn — ^ i^) in Eq- (1551) is given as 

f a + b d + e f — g c \ 
—d — e —h — i —c —f — g 
—f + g —c i — h —d + e 
V c f + g d — e a — b J 



(Al) 
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where 

2£S' {uj-r]' + r])^-{S + S'y' ^ ' 

M ^ 2EE' {uj - 1]' ^ 7]f - {E + E'Y' ^ ' 
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U — T]' + T] 




2EE' 
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-r]' + riY -{E + 


E'f 



^ M ^ 2EE' {oj - r,' + r,Y - {E + E'f ^ ' 

_ 1 ^ EE' + i'l E + E' 

2EE' {oj - i + r^y - {E + E'f ^ ' 

Here, we used the notations E = Ek, E' = Ek+q, I = Ik, I' = Ik+q, V = Vk, and r]' = r]k+q- 
For simphcity, we took to be reaL We calculate the matrix elements Eq. (IA2l) - (lA10p in 
the limit oi qd and oj / J <^ 1. Expanding Eq. f lA2p to the second order in q and w, a is 
obtained as 



A2 

AM 2^ 



1 2JEsm{kd) cosimvd) n 1 , J cosikd) smimvd) 



a ~ 

4M 

k 

, ,.nf.PAlsm'^{kd)cos'^{mvd) J'^ cos'^{kd) sm'^{mvd) 3Jf cosf/crf) cosfrnt"^,, . . , 

— - — '-^ — ^ — g5 — - + ^ — -^^r^—^U) 

By carrying out the integration over k, we obtain 

No Nqoj"^ NqJ cos{kFd) sin{mvd) 
-T-12A! + ^'" 3A^ 

, ^^2^T (3Ag - 16 J/i + V) cos(2mt;rf) , Ag - 16J^ sin^(mt;rf) 

"^^^^ ASAlcos^mvd) ^^^^ A8Alcos^{mvd) ' ^^^^^ 

Here, we have used the standard approximation, ^ j^^^j^dk F{k) -+ Nq f^^d^ F{1), where 
F{k) is an arbitrary function. By calculating Eqs. (lASp -f lATOj) in the same way as Eq. (lA2p . 
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we obtain 



6-0, (A13) 

A^'o Nqu'^ , NqJ cos(kpd) sm(mvd) 
c — TTT + qdu ^ — ^ 

^^q^^2 '^A^f^-J) + 8J'^^\ravd) _ (^^)2 (16J^ + Ag)sin^(mt;d) ^^^^^ 



A'o 

~ J ~ _ 2qdJ cos^kpd) sm{mvd)], (A15) 

, tan(mt;(i) , ,, n (2J - yu)[3 + 2tan^(mt;)] /..^x 

g ^ _^^^tanCm.d) _ (2J - ^)|3 -^2tau^(m.)l ^ ^^^^^ 

1 A'o A'ow^ , 2NoJ cosikpd) sm(mvd) 
h ^ 77 + t; TTinr + 



t/ 2 6A?, ^ 3A 



2 13A^ sin^(mt>(i) — 16 sin^(mw(i) 



24 A2 cos'^ (mvd) 
3AI + 2 J2 _ 8 + 2/^2 



-{qdYN, \' ^ , ^ , (A18) 

12A^ cos'^{mvd) 

qdtanimvd) /^,^ l\ 
I ^ ^ >- l^N, + -J . (A19) 

Note that when we calculate we have used the gap equation [Eq. fHOj) ] to 

eliminate the divergence. Assuming that the AB mode has a linear dispersion relation in 
the long- wavelength limit [qd ^ 1) and calculating the pole of Eq. fl60|) by using Eqs. 
flIl2ll - (lXT9l) . we obtain Eq. (EH]). 
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